Estimación de profundidad de disco de Secchi para el desarrollo de un algoritmo mediante técnicas de teledetección satelital

Autor
Afiliación

Mauricio Acosta

GISTAQ (UTN-FRRe)

Fecha de publicación

23 de septiembre de 2025

Resumen

Este sitio web contiene cuestiones relacionadas a la estimación de la profundidad de disco de Secchi mediante técnicas de sensado remoto. Se presenta información proveniente de investigaciones previas y se discute la confección de modelos de regresión mediante código de programación así como su evaluación mediante pruebas estadísticas.

Palabras clave

GISTAQ, UTN, FRRe, Quarto

1 Revisión bibliográfica

Esta revisión tiene como objetivo identificar los enfoques más comunes en la literatura, los tipos de datos utilizados y las variables espectrales que se consideran más relevantes para este tipo de estimaciones.

1.1 Profundidad de disco de Secchi

La profundidad de visibilidad del disco de Secchi \((SDD)\) es una medida ampliamente utilizada para evaluar la transparencia del agua. Se define como la profundidad en metros a la cual un disco blanco circular desaparece al dejarlo descender en un cuerpo de agua.

Esta variable es un indicador indirecto de la calidad del agua, ya que se relaciona con la concentración de sedimentos en suspensión, fitoplancton y otros materiales particulados. Continua siendo utilizado debido a la simplicidad y rango de aplicación universal del método, además de tratarse de un parámetro fácilmente entendible por el público.

1.2 Métodos tradicionales de estimación

En las últimas décadas, se han desarrollado múltiples ecuaciones empíricas y modelos basados en datos in situ y sensores remotos para estimar esta profundidad de manera eficiente y a gran escala. El objetivo de recopilar es identificar patrones comunes, variables predictoras relevantes y estrategias metodológicas empleadas.

Tabla 1: Algoritmos publicados para la predicción de \(\mathrm{SDD}\) utilizando la plataforma Landstat [1].
Nombre Variable Formula Muestras R²
Dekker and Peters \(\mathrm{ln(SDD)}\) \(\mathrm{ln(Red)}\) 15 0.86
Dominguez Gomez et al. \(\mathrm{SDD}\) \(\mathrm{(Green)^x}\) 16 0.90
Giardino et al. \(\mathrm{SDD}\) \(\mathrm{Blue/Green}\) 4 0.85
Kloiber et al. \(\mathrm{ln(SDD)}\) \(\mathrm{Blue/Red + Blue}\) 374 0.93
Lathrop and Lillesand \(\mathrm{ln(SDD)}\) \(\mathrm{Green}\) 9 0.98
Mancino et al. \(\mathrm{SDD}\) \(\mathrm{Red/Green + Blue/Green + Blue}\) 60 0.82

En la mayoría de los casos la relación entre \(SDD\) y la intensidad de la luz es no lineal, por lo que se utiliza \(ln(SDD)\) para realizar la regresión.

La correlación con la banda roja puede explicarse causalmente por la correlación positiva directa entre la reflectancia en el rojo y la carga bruta de partículas que induce la dispersión de partículas. De manera que mientras que la claridad del agua \((SDD)\) desciende, la intensidad en el rojo aumenta [2].

Como referencia se presentan las propiedades de las distintas bandas de la plataforma espacial Sentinel-2:

Tabla 2: Propiedades de las bandas S2-MSI, para las plataformas S2A y S2B.
Sentinel-2A Sentinel-2B
Banda Resolución espacial (m) Longitud de onda (nm) Ancho de banda (nm) Longitud de onda (nm) Ancho de banda (nm)
\(\mathrm{B01}\) (aerosol) 60 442.7 20 442.3 20
\(\mathrm{B02}\) (blue) 10 492.7 65 492.3 65
\(\mathrm{B03}\) (green) 10 559.8 35 558.9 35
\(\mathrm{B04}\) (red) 10 664.6 38 664.9 31
\(\mathrm{B05}\) (red edge) 20 794.1 14 703.8 15
\(\mathrm{B06}\) 20 748.5 14 739.1 13
\(\mathrm{B07}\) 20 782.8 19 779.7 19
\(\mathrm{B08}\) (NIR) 10 832.8 105 832.9 104
\(\mathrm{B8A}\) 20 864.7 21 864.0 21
\(\mathrm{B09}\) 60 945.1 19 943.2 20
\(\mathrm{B10}\) 60 1373.5 29 1376.9 29
\(\mathrm{B11}\) (SWIR 1) 20 1613.7 90 1616.4 94
\(\mathrm{B12}\) (SWIR 2) 20 2292.4 174 2185.7 184

1.3 Regresión lineal

En estadística, la regresión lineal es un modelo matemático usado para aproximar la relación de dependencia entre una variable dependiente \(Y\) con \(m\) variables independientes \(X_i\). Este modelo puede ser expresado como:

\[ Y = \beta_0 + \beta_1 X_1 + \cdots + \beta_m X_m\]

donde:

  • \(Y\) es la variable dependiente o variable de respuesta.
  • \(X_1,X_2,...X_m\) son las variables explicativas, independientes o regresoras.
  • \(\beta_0,\beta_1,\beta_2,...\beta_m\) son los parámetros del modelo, miden la influencia que las variables explicativas tienen sobre el regrediendo.

1.3.1 Selección de variables

La seleccion de variables y características es el foco de mucha investigación. El objetivo de la selección de variables es mejorar el rendimiento de predicción de los predictores y proporcionar una mejor comprensión del proceso subyacente que generó los datos.

La identificación de las variables más relevantes para incluir o excluir en un modelo predictivo es un paso fundamental en cualquier investigación rigurosa, especialmente cuando se busca construir modelos con alto poder explicativo y capacidad de generalización.

Comprender qué variables contribuyen en mayor medida a la calidad de las predicciones permite no solo mejorar el desempeño del modelo, sino también interpretar con mayor claridad los fenómenos que se están estudiando.

No existe una forma segura de definir la importancia de una variable puesto que la utilidad de la misma depende del modelo implementado y de las demás variables con las que interactua: una variable que es completamente inutil por si misma puede resultar en una mejora del rendimiento significativa cuando es considerada junto con otras variables [3].

A pesar de las dificultades existen estrategias útiles para determinar un subconjunto de variables útiles para la predicción. Tal es el caso de los métodos forward selection y backwards elimination. En forward selection, las variables son progresivamente incorporadas, evaluando el modelo al paso de cada una, mientras que en backwards elimination uno empieza con el conjunto entero de variables y progresivamente elimina las menos prometedoras.

Una forma de implementación de forward selection es por empezar por un modelo exento de variables independientes, pero cuyo término independiente sea la media de las variables de respuesta. A partir de este valor se calcula el residuo como la diferencia entre el valor verdadero \(Y\) y el valor predecido \(Y_{pred}\).

\[r = Y - Y_{pred}\]

Seguidamente, se computa la correlación entre cada variable independiente \(X_i\) con el residuo y se incorpora al modelo la variable con la mayor correlación absoluta. El proceso se repite esta vez con las predicciones del nuevo modelo y se continua hasta alcanzar una desempeño deseado o un número determinado de variables. Este método fue utilizado satisfactoriamente en investigaciones previas [4].

Fig. 1: Matriz de correlación entre las distintas bandas.

Una herramienta útil para ayudar en la selección de variables es una matriz de correlación. Dos variables prefectamente correlacionadas resultan redundantes en el sentido de que añadir ambas no aporta información adicional. Sin embargo, una correlación muy alta entre variables (o anti-correlación) no significa ausencia de complementariedad.

1.3.2 Generación de características lineales y no lineales

Incluir nuevas variables a partir de las originales, como productos, cocientes o transformaciones no lineales (por ejemplo, logaritmos, o raíces cuadradas) permite capturar relaciones más complejas entre las variables independientes y la variable objetivo. Por ejemplo, si se dispone de datos espectrales, pueden construirse razones entre bandas (band ratios) o índices espectrales que resalten ciertas propiedades físicas del fenómeno estudiado. Estas nuevas características permiten al modelo lineal aproximar mejor la relación entre los datos y la salida, especialmente cuando la relación real no es perfectamente lineal.

1.3.3 Métricas de desempeño

Las métricas de desempeño permiten cuantificar qué tan bien un modelo se ajusta a los datos y predice resultados. Algunas métricas evalúan la calidad del ajuste, otras penalizan la complejidad y algunas miden directamente el error de predicción.

El Coeficiente de Determinación Ajustado \(R^2-ADJ\) es una variante de R² que corrige la tendencia de este último a aumentar al añadir más variables al modelo, incluso si no aportan mejora real. Mientras que el R² tradicional mide la proporción de variabilidad de la variable dependiente explicada por el modelo, el ADJ-R²penaliza la inclusión de predictores irrelevantes, ajustando el valor según el número de observaciones y el número de variables independientes. Esto lo hace más útil para comparar modelos con distinto número de predictores, evitando la falsa impresión de mejora solo por complejidad adicional.

El Criterio de Información de Akaike (AIC) es una métrica utilizada para comparar modelos estadísticos en términos de ajuste y complejidad. Se basa en la teoría de la entropía y busca identificar el modelo que mejor explica los datos con el menor número posible de parámetros. Un valor de AIC más bajo indica un mejor equilibrio entre ajuste y simplicidad. La fórmula general es:

\[AIC = 2k - 2 ln(L)\]

donde \(k\) es el número de parámetros estimados en el modelo y \(L\) es el valor máximo de la función de verosimilitud. A diferencia del ADJ-R² el AIC no tiene un valor máximo fijo y es más útil cuando se comparan múltiples modelos candidatos sobre el mismo conjunto de datos.

El Error Cuadrático Medio (MSE) mide el promedio de los cuadrados de las diferencias entre los valores observados y los predichos por el modelo. Es una métrica sensible a errores grandes, ya que los eleva al cuadrado, y por eso penaliza más las predicciones alejadas de la realidad. Su fórmula es:

\[\text{MSE} = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2\]

donde \(y_i\) y son \(\hat{y}_i\) los valores reales y los valores predichos. Un MSE más bajo indica un mejor desempeño del modelo, pero su escala depende de la magnitud de la variable objetivo, por lo que no siempre es directamente interpretable sin contexto.

1.4 Índices espectrales

Los índices espectrales son combinaciones matemáticas de reflectancias en distintas bandas del espectro electromagnético obtenidas por sensores remotos, como satélites o drones. Estos índices realzan ciertas características o propiedades de la superficie observada, facilitando la identificación y monitoreo de elementos como vegetación, agua, suelo o áreas urbanas. Al reducir la información a una sola variable combinada, permiten detectar cambios y evaluar condiciones ambientales con mayor precisión que el análisis de bandas individuales.

1.4.1 NDWI

El Índice de Diferencia Normalizada del Agua (NDWI) es un índice espectral diseñado para resaltar la presencia y extensión de cuerpos de agua en imágenes satelitales. Se calcula utilizando reflectancias en bandas específicas del verde y el infrarrojo cercano, y su fórmula clásica es:

\[NDWI = \frac{R_{\text{verde}} - R_{\text{NIR}}}{R_{\text{verde}} + R_{\text{NIR}}}\]

donde $R_{verde} es la reflectancia en la banda verde y \(R_{NIR}\) la reflectancia en la banda del infrarrojo cercano. Valores altos de NDWI indican presencia probable de agua, ya que el agua absorbe fuertemente en el NIR y refleja más en el verde.

2 Regresión Lineal

La construcción del modelo se lleva a cabo utilizando la versión de Python 3.12 y las librerías previamente especificadas, detallando la función de cada pieza de código, observaciones, avances y descubrimientos.

2.1 Librerías

Para regresión lineal se usan principalmente las librerías de Python numpy, pandas y Scikit-learn.

  • Scikit-learn: está en el corazón de las operaciones de ciencias de datos en Python. Ofrece módulos para procesamiento de datos, aprendizaje supervisado y no supervisado, selección y validación de modelos, y métricas de error.
  • Pandas: especializada en la manipulación y el análisis de datos. Ofrece estructuras de datos y operaciones para manipular tablas numéricas y series temporales.
  • NumPy: provee al usuario con arreglos multidimensionales, junto a un gran conjunto de funciones para operar en estos.

2.2 Lectura de datos

Para la carga, lectura y manipulación de la información la librería Pandas permite convertir un archivo .csv en un DataFrame, estructura de datos que facilita la manipulación de estos. El fragmento de código siguiente agrupa los valores de \(SDD\) con los valores de reflectancia a diferentes longitudes de onda para los cuales coinciden fecha, latitud y longitud y crea un nuevo DataFrame.

import pandas as pd

# Rutas a los archivos CSV
archivo_reflectancias = "datos\\base_de_datos_gis.csv"  # contiene: fecha,punto,pixel,banda,reflect,longitud,latitud
archivo_parametros = "datos\\base_de_datos_lab.csv"     # contiene: fecha,longitud,latitud,param,valor

# Leer los archivos
df_reflect = pd.read_csv(archivo_reflectancias)
df_param = pd.read_csv(archivo_parametros)

# Filtrar los parámetros "secchi"
df_secchi = df_param[df_param["param"].str.lower() == "secchi"]

# Merge por fecha y coordenadas
merged = pd.merge(
    df_secchi,
    df_reflect,
    on=["fecha", "latitud", "longitud"],
    how="inner"
)

# Pivotear la tabla para poner bandas como columnas
tabla_final = merged.pivot_table(
    index=["param", "fecha", "longitud", "latitud", "valor"], 
    columns="banda",
    values="reflect"
).reset_index()

# Reordenar columnas: param | B01 | B02 | ... | B8A
bandas = ['B01', 'B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B11', 'B12', 'B8A']
columnas_finales = ['valor'] + bandas

# Crear Tabla final
df = tabla_final[columnas_finales]

# Guardar excel para pruebas finales
salida_excel = ".//datos//tabla_final.xlsx"
df.to_excel(salida_excel, index=False)

El DataFrame obtenido se asemeja a la siguiente tabla:

Tabla 3: DataFrame de los datos recolectados.
\(\mathrm{SDD}\) \(\mathrm{B01}\) \(\mathrm{B02}\) \(\cdots\) \(\mathrm{B12}\) \(\mathrm{B8A}\)
10 0,1728 0,1754 \(\cdots\) 0,1404 0,1869
15 0,1497 0,17022 \(\cdots\) 0,1113 0,1567
\(\vdots\) \(\vdots\) \(\vdots\) \(\cdots\) \(\vdots\) \(\vdots\)
125 0,1571 0,1563 \(\cdots\) 0,1419 0,1436
135 0,1503 0,1591 \(\cdots\) 0,1420 0,1454

2.3 Análisis de los datos

La relación no lineal entre la penetración de la luz y la profundidad de disco de Secchi fue ya descrita por diversos autores que hallan una mejor descripción de esta como una de tipo logaritmica [5]. Se verifica este comportamiento con las bandas B04, B05 y B06:

2.4 Desarrollo del modelo

Para el desarrollo del modelo se implementó la clase RegresionLineal, que contiene el código necesario para la construcción y selección de variables así como para la validación cruzada del modelo. La misma requiere unicamente de los arrays correspondientes a la variable objetivo y a las características con las que estimarla. El script se encuentra disponible acá.

Se realiza la búsqueda de ecuaciones utilizando los valores de reflectancia corregidos por SEN2COR y por ACOLITE de manera separada.

2.4.1 Corrección por SEN2COR

La implementación del método sepwise forward selection permite evaluar la adición de variables al modelo predictivo en una instancia inicial donde no se general características adicionales. Los resultados de cada paso se encuentran resumidos en la siguiente tabla:

Paso Variable R² ADJ-R² RMSE AIC
0 Ninguna -0.0263 -0.0263 44.7372 119.03
1 B05 0.2761 0.2221 37.6319 113.56
2 B02 0.7158 0.6701 23.2318 99.69
3 B08 0.7614 0.6989 21.3544 99.41
4 B03 0.8165 0.7459 18.7613 97.38
5 B04 0.8296 0.7383 17.9258 97.47

La combinación de B05, B02 y B08 resulta en el mejor ajuste para valores mínimos de RMSE y AIC.

Haciendo uso de la relación logaritmica entre la profundidad de disco y la penetración de la luz descrita y probada previamente se obtiene un modelo con una notable mejoría:

Paso Variable R² ADJ-R² RMSE AIC
0 Ninguna -0.1888 -0.1888 48.1871 121.19
1 B05 0.2949 0.2421 36.9801 112.79
2 B02 0.8333 0.8063 17.7159 91.89
3 B08 0.8379 0.7950 17.4478 93.32
4 B11 0.8519 0.7944 16.5047 93.56
5 B07 0.8549 0.7768 16.5084 95.74

La mejoría se extiende hasta la adición de la variable B07, después de la cual hay una ligera disminución de R²-ADJ, pero no lo suficientemente significativa como para tratarse de un caso de sobreajuste considerable.

La aplicación del logaritmo puede extenderse a las variables mismas, resultando en un modelo predictivo mucho más certero con una menor cantidad de variables:

Paso Variable R² ADJ-R² RMSE AIC
0 Ninguna -0.1888 -0.1888 48.1871 121.19
1 log(B05) 0.4339 0.3916 32.8131 108.64
2 log(B02) 0.8834 0.8644 15.0118 87.21
3 B06 0.8934 0.8651 14.3537 87.65
4 B08 0.9014 0.8631 13.6956 88.17
5 B12 0.9093 0.8606 13.0846 88.57

La combinación de log(B05) y B02 resulta prometedora: el valor de AIC es mínimo, el RMSE es bajo, y el R² alto. Los coeficientes y ordenada de la ecuación resultante son:

Paso Variable R² ADJ-R² RMSE AIC
0 Ninguna -0.1888 -0.1888 48.1871 121.19
1 log(B05) 0.4339 0.3916 32.8131 108.64
2 log(B02) 0.8834 0.8644 15.0118 87.21

Los coeficientes y ordenada resultantes de considerar a las primeras dos variables son:

Variable Coeficiente
log(B05) -12.260673
log(B02) 13.030511
Ordenada 5.275337

El uso de cocientes entre las variables otorga una cantidad mucho más extensiva de predictores, observamos que de considerar a estos junto con las variables linealizadas, estos predominan como mejores contribuyentes a la mejora de la predicción:

Paso Variable R² ADJ-R² RMSE AIC
0 Ninguna -0.1888 -0.1888 48.1871 121.19
1 B05/B02 0.8822 0.8734 15.0707 84.77
2 log(B01) 0.8966 0.8798 14.0689 84.36
3 B08/B06 0.9091 0.8850 13.1732 84.66
4 B8A/B11 0.9100 0.8750 13.0816 86.29
5 B07/B11 0.9087 0.8596 13.1518 88.25

Los coeficientes y ordenada resultantes de considerar a las primeras dos variables son:

Variable Coeficiente
B05/B02 -4.268558
log(B01) 0.643369
Ordenada 9.430494

2.4.2 Corrección por ACOLITE

Se repite el proceso que con los datos corregidos por el algoritmo SEN2COR. Se menciona que en la lectura y tratamiento de datos se descartaron los puntos para los cuales los valores de reflectancia eran negativos luego de la corrección.

Partimos de haber hallado que la aplicación del logaritmo a la profundidad de disco resulta en un mejor ajuste. En relación a los valores corregidos por ACOLITE esto devuelve los siguientes resultados:

Paso Variable R² ADJ-R² RMSE AIC
0 Ninguna -0.2303 -0.2303 47.8282 59.46
1 B05 0.1189 -0.0378 38.8925 56.79
2 B02 0.6160 0.4521 27.3299 51.54
3 B08 0.7201 0.4925 23.1573 51.70
4 B03 0.8162 0.5245 16.6068 48.02
5 B06 0.8334 0.2594 16.2201 49.84

Observamos que para las mismas combinaciones los valores de R² y R²-ADJ aumentaron significativamente respecto a aquellos con la corrección por SEN2COR.

Para una relación doblemente logarítimica, obtenemos lo siguiente:

Paso Variable R² ADJ-R² RMSE AIC
0 Ninguna -0.2303 -0.2303 47.8282 59.46
1 log(B05) 0.3287 0.2083 28.9200 52.62
2 log(B02) 0.6888 0.5550 21.0484 49.21
3 B01 0.7729 0.5845 18.4856 48.91
4 log(B08) 0.8185 0.5269 15.6606 48.45
5 B06 0.8315 0.2347 15.0506 50.27
Variable Coeficiente
0 log(B05) -5.058134
1 log(B02) 6.861923
2 B01 -18.872053
3 Ordenada 7.577205

Donde esta vez se ve una desmejora respecto al modelo corregido por SEN2COR para un modelo de dos variables logarítimicas.

La incorporación de cocientes al abanico de variables devuelve las siguientes métricas:

Paso Variable R² ADJ-R² RMSE AIC
0 Ninguna -0.2303 -0.2303 47.8282 59.46
1 B04/B02 0.7325 0.6853 20.3765 47.21
2 B04/B01 0.7795 0.6848 19.1909 48.01
3 B07/B06 0.8216 0.6737 15.8882 47.67
4 B08 0.8167 0.5294 16.0595 49.33
5 B11/B05 0.8592 0.3978 13.5437 48.60
Variable Coeficiente
0 B04/B02 -1.582140
1 Ordenada 6.402011

Nuevamente, se ve una desmejora respecto al modelo corregido por SEN2COR. El número de variables se reduce a una.

2.5 Ecuaciones halladas

El par de mejores ecuaciones halladas para cada algoritmo de corrección utilizando la transformación logarítmica fueron las siguientes: \[ log(SDD)_{SEN2COR} = -12.26 \cdot log(B05) + 13.03 \cdot log(B02) + 5.27 \]

\[ log(SDD)_{ACOLITE} = -5.06 \cdot log(B05) + 6.86 \cdot log(B02) -18.87 \cdot B01 \]

Sus ajustes se visualizan de la siguiente manera:

El par de mejores ecuaciones halladas para cada algoritmo utilizando cocientes entre las variables es:

\[ log(SDD)_{SEN2COR} = -4.27\frac{B05}{B02} + 0.64\cdot log(B01)+ 9.43 \]

\[ log(SDD)_{ACOLITE} = -1.58\frac{B04}{B02} + 6.40 \]

Sus ajustes se visualizan de la siguiente manera:

Observamos que la utilización del algoritmo de ACOLITE para la corrección atmosférica ocasiona una desmejora en la estimación.

2.6 Mapas de la estimación

Los mapas generados para cada fecha corresponden al recorte del área de interés del producto satelital corregido por ACOLITE.

3 Métodos de aprendizaje automático

[1]
H. Rubin et al., «Remote Sensing of Lake Water Clarity: Performance and Transferability of Both Historical Algorithms and Machine Learning», Remote Sensing, vol. 13, p. 1434, abr. 2021, doi: 10.3390/rs13081434.
[2]
M. Matthews, «A current review of empirical procedures of remote sensing in Inland and near-coastal transitional waters», International Journal of Remote Sensing - INT J REMOTE SENS, vol. 32, pp. 1-45, nov. 2011, doi: 10.1080/01431161.2010.512947.
[3]
I. Guyon y A. Elisseeff, «An Introduction to Variable and Feature Selection», Journal of Machine Learning Research, vol. 3, pp. 1157-1182, 2003, Disponible en: https://www.jmlr.org/papers/volume3/guyon03a/guyon03a.pdf
[4]
M. Bonansea et al., «Evaluating the feasibility of using Sentinel-2 imagery for water clarity assessment in a reservoir», Journal of South American Earth Sciences, vol. 95, p. 102265, 2019, doi: https://doi.org/10.1016/j.jsames.2019.102265.
[5]
G. Wu, J. De Leeuw, A. Skidmore, H. Prins, y Y. Liu, «Comparison of MODIS and Landsat TM5 images for mapping tempo-spatial dynamics of Secchi disk depths in Poyang Lake National Nature Reserve, China», International Journal of Remote Sensing 29 (2008) 8, vol. 29, abr. 2008, doi: 10.1080/01431160701422254.